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Abstract 

We introduce an age-structured asexual population model containing 
all the relevant features of evolutionary ageing theories. Beneficial as 
well as deleterious mutations, heredity and arbitrary fecundity are present 
and managed by natural selection. An exact solution without ageing is 
found. We show that fertility is associated with generalized forms of the 
Fibonacci sequence, while mutations and natural selection are merged into 
an integral equation which is solved by Fourier series. Average survival 
probabilities and Malthusian growth exponents are calculated indicating 
that the system may exhibit mutational meltdown. The relevance of the 
model in the context of fissile reproduction groups as many protozoa and 
coelenterates is discussed. 
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1 Introduction 



For all alive organisms death is inexorable. If life is not abbreviated by 
disease, predation or accidents then death is preceded by senescence ||. 
Senescence seems to be related to reproductive strategies and it starts after 
an individual reaches the fertile age. It can be a slow process for species which 
reproduce many times (iteroparous) or an abrupt process for species reproducing 
only once (semelparous) - the catastrophic senescence of the Pacific salmon being 
a good example of the late. 

There arc tree kinds of ageing theories: biochemical, evolutionary ^, ^ 
and telomeric The biochemical invokes damages on DNA, cells, tissues 
and organs. Defective proteins are synthesized altering the normal course of 
metabolism. The presence of free radicals, i. e., of unpaired highly reactive elec- 
trons, can cause death of the cells or may even lead to cancer. As a consequence, 
such theories are today firmly connected to modern gerontology. Biochemical 
theories predict that species with higher metabolism would have shorter lifetime. 
A criticism against this result arises when the lifetime of birds and mammals (of 
the same size and under optimal conditions in captivity) are compared. Usually, 
birds live longer. Even closely related organisms like bats and rodents with com- 
parable sizes have very different lifetimes. Such differences should be explained 
by the evolutionary theories of ageing. In the telomere hypothesis of senescence, 
replication of normal cells is accompanied by a telomeric shortening. This acts 
as a mitotic clock resulting in a permanent exit of the cell cycle. 

Evolutionary theories of ageing are hypothetico-deductive in character, not 
inductive. They fall into two classes: the optimality theory and the mutational 
theory. In the optimality theory fitness is maximized by increasing the 

survival and reproduction rate early in life at the expense of late, i. e., it sees 
senescence as a necessary cost of processes beneficial to youth. On the other 
hand, the mutational theory |^ explains senescence as a result of late-acting 
deleterious mutations, i. e., that a greater mutation load on the last part of the 
life is less strongly selected. 

In this paper, we obtain both analytical and numerical solutions for a simple 
age-structured population model. In this model, reproduction is asexual and 
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the individuals are submitted to helpful or deleterious hereditary mutations. 
The number of offsprings is fixed but arbitrary otherwise and the population 
dynamics is managed by natural selection. Amazingly, this system does not 
exhibit senescence even when deleterious mutations are predominant. From a 
mathematical point of view, the solution wc found factorizes into a fertility 
and a mutational sectors. It is shown that the fertility sector is completely 
described by generalized forms of Fibonacci sequences. On the other hand, the 
mutational sector, solved by Fourier series, incorporates the combined effects of 
mutation and natural selection. If harmful mutation is intense then a mutational 
meltdown ^ process can take place. All these results were corroborated by 
Monte Carlo simulations. From a biological point of view, our results make the 
model a good candidate to describe groups in which all reproduction occurs by 
fission, such as protozoa and a miscellany of coelenterates, since all them appear 
to lack ageing 

2 The Model 

Consider a population distributed by L + 1 ages i (i — 0,1, L) with re- 
spective birth rates m^. We call babies the individuals at age 0. They do not 
reproduce, i. e., toq — 0. Individuals with age i = L die after reproduction. 
Let Ni{J,t) be the number of individuals at age i with survival probability J 
between J and J + dJ dX time t (of course, J G [0, 1]). We choose, as initial 
condition, an uniform distribution of babies, i. e., Ni{J,Q) — Nq Si^ with Nq 
being a constant. Let us point out that senescence here means that the average 
survival probability drops with age i. At time t, each individual is submitted 
to mutation which changes its survival probability from J to J' in the following 
way 

J' = J (1) 

where e is a random mmibcr chosen in the interval [—a,b] (a > and 6 > 0), 
generated by the uniform probability density distribution q{e) = 
{H{x) stands for the Heaviside function). If a > then deleterious mutations 
are dominant. The transition probability P{J',J), that an individual has its 
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survival probability changed from J to J' is given by the integral P{J',J) = 
J-^ (5( J' — Je'^)q{e)dt. Observing that a helpful mutation is allowed as long as 
it does not increase the survival probability beyond unity, we get 



J' {a + b)P{J', J) = {[H{J' - Je-'') ~ H{J' - Je'')][H{J) - H{J - e"^)] 

+ [H{J' - Je-") - H{J' - 1)][H{J - e-'') - H{J - 1)]}(2) 

After mutation, all population passes through natural selection. The sur- 
vivors then reproduce, generating babies with inherited characteristics, i. e., 



No{J, i) = ^ TO, iV,( J, t), fort>l 



(3) 



At time step i + 1 individuals get older or die. Taking into account mutations 
and natural selection (in this order), their number is given by 



N,+i{f,t+l)= J' P{f,J)N,{J,t)dJ fori = 0,...,L-l (4) 
Jo 

Clearly, Ni{J, i < i) = for « ^ 0. For t > 1, equations (3) and (4) can be 
rewritten in the matricial form 



N{J', t+l) = M / /(J', J) N{J, t) dJ 
Jo 



(5) 



where f{J',J) — J' P{J',J) and N{J,t) and M are the following vector and 
matrix 



NiJ,t) = 



( N,{J,t) \ 

N2{J,t) 



, M = 



f nil 1712 
1 
1 



TOL-1 rriL \ 
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1 



/ 



Matrix M is a special form of the Leslie matrices Iterating equation 

(5) we obtain at time t > 1 
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where 




(7) 



and M° is the L x L identity matrix. 

3 Generalized Fibonacci Sequences 

Equation (6) shows that the dynamics factorizes into two sectors: the fertil- 
ity, exclusively contained in the matrix M, and the mutational, enclosing both 
mutations and natural selection processes and represented by the function F. 
Clearly, at time t, we need only to know the first column elements of the (i — 
power of the fertility matrix. As we shall see, these elements form a Fibonacci 
sequence. Let us call AL{t + 1 — i) the element of the ith line and first column. 
It is a simple exercise to verify that if L = 2 and mi = = 1 then A2 (t) can 
be calculated, at any time t, through the relation A2{t) ^ A2{t — 1) + A2{t — 2), 
with ^2(1) = ^2(2) = 1, which is exactly the Fibonacci sequence!. The ratio of 
two sucessives numbers, limt^oo ^L(i + l)ML(i) = 1-618 . . . , gives the golden 
section or the divine proportion as called by Kepler. It is astonishing to find 
this ubiquitous sequence in so disparate things like the cluster-cluster aggregates 
[ pi] , the division of a line into extreme and mean ratio, the pentagram star worn 
by the Pythagoreans, the continued fractions, the aesthetic proportions of the 
Parthenon at Athens ]T^ , and now, here, in population dynamics. 

When the number of ages is 4 (L = 3) and mi = m2 = = 1 then 



Asit) = Asit - 1) + A-iit - 2) + Asit - 3), with ^3(1) = ^3(2) = 1 and 



^3(3) = 2 (this sequence generates the so called tribonacci numbers). For an 
arbitrary number of ages and fecundity we have 
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L 

^iW-^mfc AL(i-fc) , fort> (i + 1)) (8) 

fe=i 

The first L numbers (necessary to initialize the sequence above) are evaluated 
through the identification Ai^{t) = A2{t), for t = 1, . . . ,L. The numbers A2{t) 
are determined, on the other hand, by using ^2(1) = 1 and ^2(2) = mi in the 
expression above. Recurrence formulas like equation (8) are generalized forms 
of the Fibonacci sequence. These results, together with equation (6), allow us 
to write down the number of individuals at time t with age i > 1 and survival 
probabibity J (we renamed Jt J) 



N,{J,t)^ Nq AL{t + l-i) F{J,t), iit>i 

= 0, otherwise (9) 

The number of babies NQ{J,t) may be determined using equations (3) and 
(9). The mean survival probability at any time t and at arbitrary age i can be 
calculated as 

N,iJ,t) dJ J^FiJ,t)dJ 
which is independent of i, that is, the model does not show senescence !. It is 
important to note that, although ageing is absent in our model (in the sense 
that the mean survival probability does not depend on the age i), an individual 
can show senile decay (that is, its survival probability J diminishes with time). 
This individual handicap is compensated by the natural selection of the fittest 
so, at a collective level, no senescence is observable. We also call attention to the 
fact that, in our model, the total number of individuals with age i diminishes 
with i. If we integrate Ni{J,t) over J (equation (9)), the ratio between two 
successives ages i and i + 1 is given by A^it + 1 ~ i) / A^it + i) of the Fibonacci 
sequence. This ratio depends on age i only for small time t. Moreover, the 
logarithm of this ratio corresponds to the so called mortality rate. For humans, 
the Gompertz law suggests that mortality increases exponentially with age. 
In our model it is constant. 
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At a fixed age i and survival probability J, the population increases with 
time at a rate given by ]^l^jft'!i\) — ^^Lit-'i)^ ' F0,t'-i) • "^^^ ^^^^ ratio of the 
right side can be easily calculated by the generalized Fibonacci sequences, but 
the second requires a numerical analysis. 

4 Numerical Analysis 

The continous variable J can be divided into Q intervals such that J = 
j /Q for j ~ 1 . . .Q and infinitesimal increment dJ = 1/Q. For Q big enough 
we expect to reobtain the continous limit. In the same way, the products of 
f{Jj,Jj-i) in the equation (7) can be seen as simple matricial products. We 
wrote down a FORTRAN program with extended precision to determine F{J,t) 
and < J >i (t) at any elapsed time t. For Q = 400, L = 10, a = 0.04 and 
b = 0.02, Fig.l shows the dependence of F{J,t) with J at different instants. 
After a time t > 50 we verified that 

F{J,t+l) = cF{J,t) (11) 

where c depends on a and b but not on J. This means that after enough time, 
the system reaches an asymptotic limit where F{J,t) = c *F{J), i. e., there is 
a separation of variables and F{J) is a stationary solution. Fixing h — 0.02 and 
varying a = 0.02, 0.04 and 0.08, we determine numerically that c = 0.92, 0.76 
and 0.48, respectively. If, for example, we choose L — \Q and nii — 1 for any i 
then limj^oo '^Ai.ti)' ~ ^-^^^^ ■ ■ ■' ^^'^ ^ Malthusian exponential growth of the 
population e*" * with r = 0.61,0.42 and —0.04, respectively, is obtained. The 
last value shows that the system exhibits mutational meltdown (extinction) . 

The discretized form of the variable J can also be used in order to calcu- 
late the mean survival probability < J{t) > (equation (10), dropping the now 
unnecessary lower index i) as a function of time. Table 1 shows our results for 
L = 10. For fixed matrix dimension, the time convergence is very fast. 

5 Monte Carlo Simulation and Exact Solution 

To check out all these features, we also made some simulations of the model. 
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Starting with an uniform distribution of the babies, we submit all of them to 
mutations as described in equation (1). For each baby, a random number e 
is generated in the computer and its new survival probability J' is calculated. 
Then, playing the role of the natural selection, a random number r is generated 
and compared with J'. If r < J' then the baby becomes an individual of age 1 
and produces mi offsprings with inherited characteristics (that is, with the same 
survival probability J'). As the process continues, care should be taken in order 
to avoid an explosion of the computer's memory. To this end, we limited the 
number of individuals by a random decimation. Usually, in population dynamics 
simulations, this decimation is interpreted as a result of food restrictions |Q. 
Fig. 2 (a) shows that for L = 10, a = 0.04 and b = 0.02, the final mean 
survival probability < J >— limt^oo < J > (t) approaches 0.78 in a very 
good agreement with our numerical results. For a = 0.02 and b — 0.02 we find 
< J >~ 0.96. The case a — 0.08 (Fig. 2 (c)) leads to extinction, as it was 
predicted. 

As another important check, let us look at the Euler-Lotka equation |^ 
which, for our model, reads 

L 

Y^iii, {<J> e-^Y (12) 

i=l 

with r being the Malthusian growth exponent. 

Substituting the simulated values < J 0.78 and < J >= 0.96 into the 
Euler-Lotka equation, we obtain the Malthusian growth exponents r = 0.44 and 
r = 0.65, respectively, which are in a fairly good agreement with those of our 
numerical prediction. 

Incidentally, equations (11) and (7) can be consistently used to guide us 
to an analytical solution of the stationary F(J) . One can write the integral 
equation 

cF(j')- / f{J\J)F{J)dJ (13) 

Integrating the right side and expanding the result in a Fourier series, equa- 
tion (13) turns out to be a set of linear equations for the Fourier coefficients. 
In order to have non trivial solutions, this set must have null determinant. The 
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condition of zero determinant allow us to obtain the constant c. For example, 
if a — 0.04 and b ~ 0.02 we find c — 0.75, in well accordance with the numer- 
ical results. The Fourier coefficients have awful expressions so, to preserve the 
reader, we will not give them. 

6 Discussion 

Let us discuss the relevance of our model by comparing it with others evolu- 
tionary models. The Pcnna model is certainly the most intensively investigated 

. It exhibits ageing and sometimes catastrofic senescence. Contrary to what 
happens in our case, in the Penna model only babies are affected by mutations. 
Moreover, mutation plays the role of a programmed death - individuals which 
have accumulated a number of mutations (i. e., number of I's in the bit-string) 
larger than a threshold T die. This fate can be anticipated only if the individual 
dies by food restrictions (Verhulst factor, which acts irrespective of individual 
fitness). In our model there is not such a threshold and individuals may live 
longer. Besides that, natural selection here operates in a very hard (and ex- 
plicit) way to eliminate those individuals which have suffered bad mutations. 
Also exact solutions have been found for the Penna model |T^. The evolution 
equations involve directly the Verhulst factor. We have similar equations ((3) 
and (4)), but with the survival probability J instead. Amazingly, the Leslie 
matrices were also found in the Penna model but in a different context. There, 
the elements of these matrices are connected to the mutation rate while in our 
model they are associated to the birth rate. 

As we said in the beginning, there are two kinds of evolutionary theories: op- 
timal and mutational. Our model belongs to the latter. The optimality theory 
is based on the fact that some genes have antagonistics effects, that is, they can 
be very beneficial early in life but deleterious late. For example, genes enhanc- 
ing early survival by promoting a bone hardening might reduce later survival 
by promoting arterial hardening. These ideas were completely embodied by 
the Partridge-Barton model 0. Further studies on this model, have incorpo- 
rated somatic as well as hereditary mutations leaving the model with 
two mechanisms of senescence: antagonistic pleiotropy and accumulation of bad 
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mutations. But ageing (due to mutations) emerges in these works as a result 
of turning more intense with age (in some artificial or arbitrary way) the muta- 
tional strength. Their procedures would be equivalent to assume, in our model, 
that the e (the mutational control parameter of eq. (1)) is a function of the age 
i. Clearly, this would also trigger an ageing process in our model. 

More interesting, however, is a different version of the Partridge-Barton 
model which is called VoUmar-Dasgupta ||l^ model. It was generalized by 
Heumann and Hotzel |^o| to support many age intervals. In this model, each 
individual carries, like its genome, a whole set of independent survival probabil- 
ities {Jo, ■ ■ ■ , Jl} where Ji (i = . . . L) is the actual survival probability at age 
i. Deleterious mutations can now affect any of them but only those coincident 
with the actual age i will pass through natural selection. This means that most 
of the damage will only manifest later on. This accumulation of harmful mu- 
tations leads to senescence. Our model differs from Heumann and Hotzel only 
in the point that our individuals carry just one survival probability - that of its 
actual age. This is sufficient to change radically the results. 

In summary, although containing all the relevant features of evolutionary 
systems like age-structure, advantageous or deleterious mutations, reproduc- 
tion with inherited characteristics and natural selection, our model does not 
show senescence. In this way, it is a good candidate to a appropiately describe 
some coelenterates and prokaryotes groups, since all them appear to lack senes- 
cence. On the other hand, the beautiful analytical solution that we find and the 
techniques involved, encourage us to look forward new and more sophisticated 
models. 
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TABLE CAPTION 



Table 1. The mean survival probability < J >, obtained by using our 
matricial formalism, as a function of the matrix dimension Q and time t for 
L = 10, a = 0.04 and b = 0.02. 
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FIGURE CAPTION 



Figure 1. Plots of the function F{J,t) for t = 100, 200, 400 and 800 against 
the survival probabihty J. In order to bring them to the same scale, they were 
multiplied by a factor s = 10^\ 10^^, 10^'^ and 10^^, respectively. 



Figure 2. (a) time evolution of the simulated values of < J > (t) for L = 10, 
a = 0.04 and b = 0.02. The eleven age curves coincide and are undistinguish- 
ablc, (b) the corresponding population plotted against the time. The stationary 
behavior is an artefact of the decimation process, (c) the same for L = 10, 
a = 0.08 and h = 0.02, (d) the corresponding population comes to extinction. 
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